function matpath =  extractMetrics(matpath,time,channels,positions,maxObjects,options)

%Thomas Stevens - Lim Lab
%opened September 2012
%-------------------------------------------------------
%matpath =  extractMetrics(matpath,time,channels,positions,maxObjects,options)
%Inputs: (Optional) matfile path with image stack, binaries, metadata structure
%specific positions, channels, or timepoints; parallelization option
%Outputs: cell array of structures {c,p} -- need cell since number of
%objects in a stack is variable.
%-------------------------------------------------------
%standard shape and pixel intensity properties.
%Custom fields calculated with previous BW masks.
%Centroid will be time-averaged. Box bounds all time.
%Future: extend with 4D watershed.

%% INPUT Defaults
if nargin==0 || ~ischar(matpath)
    [file,path] = uigetfile();
    matpath = [path file];
end
if nargin < 2 || isempty(time)
    m = matfile(matpath);
    s = m.s;
    time = 1:s.time; %all
end
if nargin < 3 || isempty(channels)
    %fluorescent only
    m = matfile(matpath);
    s = m.s;
    channels = find(ismember(s.channel_name_std,'RGB'))';
end
if nargin < 4 || isempty(positions)
    m = matfile(matpath);
    s = m.s;
    positions = 1:s.positions; %all
end
if nargin < 5 || isempty(maxObjects)
    maxObjects = 20;
end
if nargin < 6 || isempty(options)
    options = {''};
end

if ischar(options), options = {options}; end

isParallel = startParallel(options);

%Initialize matfile
m = matfile(matpath,'Writable',true);
s = m.s;

%h_bar = waitbar(0,'Quantifying Features...');

%Only nonredundant N-D properties
prop_default_cell = {'Area','Centroid','BoundingBox','Image','PixelValues', ...
    'WeightedCentroid','MeanIntensity','MinIntensity','MaxIntensity', ...
    'SubarrayIdx'};
prop_custom_cell = {'Area_t','Channel_image','Channel_area','Channel_area_t'};

%Filtering values, 3d
open_size = max(time)*ceil(20/prod([s.calibX s.calibY s.calibZ]));

%Initialize props structure in matfile
disp('Initializing property structure in matfile')
tic
props_unit = struct();
for field = cat(2,prop_default_cell,prop_custom_cell)
    field_c = char(field);
    if strcmp(field,'SubarrayIdx')
        props_unit = setfield(props_unit,field_c,cell(1));
    else
        props_unit = setfield(props_unit,field_c,[]);
    end
end
props_unit = orderfields(props_unit);
%Need to make last dimension non-singleton when pre-allocating
props_rep = repmat(props_unit,[maxObjects max(channels) max([positions 2])]);
m.props = props_rep;
%clear props_unit -- keep for null indexing
toc

%Needs to be maximal except on singleton dimensions, where all 0
conn = conndef(ndims(cell2mat(m.BW(2,1))),'maximal');
%Component Connectivity: No z-neighborhood if 2d in time (yields border errors)
if s.z==1
    conn(:,:,1,:)=0;
    conn(:,:,3,:)=0;
end

%Moved to read_nd2()
%Transient-clearing connectivity: minimal time only
%time_border_conn = zeros(3,3,3,3);
%time_border_conn(2,2,2,:)=1;


tic
tempDir = [matpath(1:find(matpath=='/',1,'last'))  'workers/'];
mkdir(tempDir)
baseFile = [tempDir 'part_%03d.mat'];
parfor p=positions
    %for p=positions
    disp(['Clearing border and transients, Position #' num2str(p)])
    tic
    %Re-point to matfile for parallelization
    m = matfile(matpath,'Writable',true);
    if isParallel
        %Workers save to distinct files in tempdir
        m_part = matfile(sprintf(baseFile, p),'Writable',true);
        m_part.props = props_rep;
    else
        %One core read/write: save directly to mat
        m_part = matfile(matpath,'Writable',true);
    end
    for c=channels
        if s.channel_name_std(c)=='R'
            BW_cp = cell2mat(m.BW(c,p));
            BW_p_r = BW_cp(:,:,:,time);
        elseif s.channel_name_std(c)=='G'
            BW_cp = cell2mat(m.BW(c,p));
            BW_p_g = BW_cp(:,:,:,time);
        else
            continue
        end
    end
    toc
    disp('Identifying XYZT-connected components')
    tic
    %Cluster masks from combined Red/Green fluorescent channels -- OK
    %Need to include, validate watershedding
    CC = bwconncomp(BW_p_r | BW_p_g, conn);
    %% Filtering
    %clear BW_cp BW_p_r BW_p_g
    [BW_cp,BW_p_r,BW_p_g] = deal([],[],[]);
    L = labelmatrix(CC);
    filter_props = regionprops(L,'Area');
    idx = find([filter_props.Area] > open_size);
    BW_filt = ismember(L,idx);
    %Consider thicken Inf --> and with next to eliminate joinings
    %clear L CC filter_props idx
    [L,CC,filter_props,idx] = deal([],[],[],[]);
    %Load a position of raw image data
    I_p = cell2mat(m.I(p,1));
    for c=channels
        disp(['     Computing properties for C=' num2str(c) '...'])
        %get all N-D builtin region properties
        %Need to index on time for every time slice
        props_cp = regionprops(BW_filt,I_p(:,:,:,time,c),prop_default_cell);
        %% Add custom properties
        %Area_t = timecourse of total image Area
        %Channel_image = Image AND BW{c,p}
        %Channel_area = Area (all dims) of Channel_image
        %Channel_area_t = timecourse of Channel_area
        BW_cp = cell2mat(m.BW(c,p));
        BW_cp = BW_cp(:,:,:,time);
        %for i=1:length(props_cp) -- throws error
        for i=1:min(maxObjects,length(props_cp))
            props_icp = props_cp(i);
            image_i = props_icp.Image;
            %Sum binary pixels in bounded image over all x,y,z
            props_icp.Area_t = squeeze(sum(sum(sum(image_i,1),2),3));
            %Take only part of mask in selected conncomp
            idx = props_icp.SubarrayIdx;
            c_image_i = image_i & BW_cp(idx{:});
            props_icp.Channel_image = c_image_i;
            %props{c,p}(i).Channel_raw = image_i & Icp(idx{:});
            %Sum binary pixels in bounded image over all x,y,z for one
            %channel
            props_icp.Channel_area = sum(c_image_i(:));
            props_icp.Channel_area_t = squeeze(sum(sum(sum(c_image_i,1),2),3));
            %write to structure array in matfile
            %disp(['i=' num2str(i) ',c=' num2str(c) ',p=' num2str(p)])
            m_part.props(i,c,p) = orderfields(props_icp);
        end
        %Status
        %if ~isParallel
        %   waitbar( sum([p-1 c].*[channels 1])/prod([channels positions]), h_bar)
        %end
    end
    toc
end
toc

if isParallel
    %cat separate worker files
    tic
    disp('Loading temp worker files into central file')
    m = matfile(matpath,'Writable',true);
    for p=positions
        m_part = matfile(sprintf(baseFile, p));
        for c=channels
            disp(['Position #' num2str(p) ', Channel #' num2str(c)])
            for i=1:maxObjects
                m.props(i,c,p) = m_part.props(i,c,p);
            end
        end
    end
    %save s structure
    m.s = s;
    rmdir(tempDir,'s')
    toc
else
    %close(h_bar)
end

end